Process-oriented analysis of dominant sources of uncertainty in the land carbon sink

The observed global net land carbon sink is captured by current land models. All models agree that atmospheric CO2 and nitrogen deposition driven gains in carbon stocks are partially offset by climate and land-use and land-cover change (LULCC) losses. However, there is a lack of consensus in the partitioning of the sink between vegetation and soil, where models do not even agree on the direction of change in carbon stocks over the past 60 years. This uncertainty is driven by plant productivity, allocation, and turnover response to atmospheric CO2 (and to a smaller extent to LULCC), and the response of soil to LULCC (and to a lesser extent climate). Overall, differences in turnover explain ~70% of model spread in both vegetation and soil carbon changes. Further analysis of internal plant and soil (individual pools) cycling is needed to reduce uncertainty in the controlling processes behind the global land carbon sink.

The observed global net land carbon sink is captured by current land models. All models agree that atmospheric CO 2 and nitrogen deposition driven gains in carbon stocks are partially offset by climate and land-use and land-cover change (LULCC) losses. However, there is a lack of consensus in the partitioning of the sink between vegetation and soil, where models do not even agree on the direction of change in carbon stocks over the past 60 years. This uncertainty is driven by plant productivity, allocation, and turnover response to atmospheric CO 2 (and to a smaller extent to LULCC), and the response of soil to LULCC (and to a lesser extent climate). Overall, differences in turnover explain~70% of model spread in both vegetation and soil carbon changes. Further analysis of internal plant and soil (individual pools) cycling is needed to reduce uncertainty in the controlling processes behind the global land carbon sink.
Over the last 60 years, there has been a continuous rise in anthropogenic CO 2 emissions. Around the equivalent of a quarter of these emissions have been taken up by the land biosphere (known as the natural land sink), acting as a strong negative feedback to mitigate climate change 1 . To be able to project the global carbon cycle (and hence the climate response) in the future, we need to understand the underlying processes and their timescales that drive the contemporary land sink. There are many distinct but interdependent mechanisms that regulate the flow of carbon into, through, and out of the land. The timescales at which these processes (e.g. photosynthesis, allocation, plant growth, litterfall, plant mortality and soil turnover) act on the carbon cycle range from days to centuries, and the interplay between these determines the changes in land carbon storage 2 .
Recent work suggests the global net land sink is located in northern latitudes 3 . The attribution of drivers is highly uncertain, but a combination of increasing atmospheric CO 2 concentrations 4 , reactive nitrogen (N) deposition 5 and atmospheric warming 6 are likely responsible. Tropical lands are probably nearer net zero carbon sinks due to large land-use and land cover changes (LULCC) carbon losses counteracting the 'natural' sink 7,8 . Atmospheric observations also show that the Amazon forest is currently carbon neutral ('natural' sink = LULCC-driven source) 9 , and this provides further evidence of a limited net tropical sink. Yet, the largest gross fluxes between the land and atmosphere are often found in tropical regions 10 , and so changes in tropical ecosystem functioning can have significant global impacts.
Process-based dynamic global vegetation models (DGVMs) that simulate processes of carbon uptake and release can help to elucidate the roles of individual drivers (rising atmospheric CO 2 , changes in climate, nutrient deposition, and LULCC), attribute the change to processes, and quantify regional sinks over timescales needed given the multitude of timescales over which terrestrial processes act (i.e. beyond the period of reliable empirical and remote-sensed data).
DGVMs are used to estimate the natural land sink (hereafter simply land sink) as part of the global carbon budget (GCB 8 ). Globally, the DGVM multi-model mean estimate of the global land sink is consistent with the global carbon budget residual land sink (the difference between fossil and land-use emissions and atmospheric and ocean sinks, see Table 5 in ref. 8), however, there is a significant spread across models, and here we show the spread widens at regional scales, when quantifying changes in vegetation and soil carbon, or attributing changes in internal processes to external drivers. Therefore, while model ensembles are helpful in analysing global-scale processes, they must be interpreted in the context of their process representation and unique and collective biases 2,4 .
These deficiencies are of high significance when trying to understand past changes but also when DGVMs are used to make predictions of future carbon cycling 11 , as limited trust can be placed in projections when certain fundamental processes and carbon-climate feedbacks are not fully captured 2 . Therefore, it is of high priority to identify the leading sources of uncertainty in DGVMs and understand the relationship between drivers and processes on varying spatial scales. In this study, we use a triple (3D-matrix) approach to identify where models agree and, just as importantly, disagree, and thus guide future modelling efforts: • What are the (1) external drivers (concurrent rises in atmospheric CO 2 and N deposition, climate and LULCC), (2) main regions (tropics and extra tropics) and (3) processes (production vs turnover) primarily responsible for the changes in the global net land carbon sink?
We use the suite of 18 DGVMs from the GCB2021 (TRENDYv10; ref. 8) to quantify changes in net carbon exchange and carbon stocks over the period 1959-2020. TRENDYv10 provides a set of simulations to attribute these changes to drivers, and we use a process attribution framework to decompose changes in carbon stocks into those driven by productivity and turnover separately. This framework enables us to express productivity and turnover-induced changes in carbon stocks in units of PgC, which allows for a direct comparison between both processes, which have units of carbon per unit time and time, respectively (see Methods).

Drivers of global and regional land sinks
The global net land sink is derived from the difference between the fossil fuel emissions (E FOS ) and the CO 2 accumulation in the atmosphere (G ATM ) and uptake by the oceans (S OCEAN ). We refer to this estimate of the net land sink as the 'observed' GCB budget constraint (= E FOS -G ATM -S OCEAN ; see methods and ref. 8) which grew from 0.2 ± 0.4 (mean ± std. dev) PgC yr −1 in the 1960s to 1.7 ± 0.6 PgC yr −1 in the decade 2011-2020, with DGVMs (S3 simulation; see Methods) capturing the increase (−0.1 ± 0.6 to 1.6 ± 0.5 PgC yr −1 ; Fig. 1). DGVMs suggest the enhanced net sink over the past 60 years has mainly been driven by rising atmospheric CO 2 concentrations and nitrogen deposition (1.2 ± 0.2 PgC yr-1 over 1960-1969 rising to 3.5 ± 0.8 PgC yr-1 over 2011-2020) (Fig. 1b), with the sink partially offset by relatively constant net LULCC emissions of 1.3 ± 0.5 PgC yr −1 (predominantly arising from tropical deforestation and shifting cultivation 12 ; Fig. 2 and Supplementary Fig. 1). DGVMs indicate that climatic variability drives the large year-to-year changes (±2 PgC yr −1 ) in the land sink ( Fig. 1), with long-term (multi-decadal) climate trends reducing the land sink since the 1980s (−0.4 ± 0.5 PgC yr −1 over 1980-2020) (Fig. 1b).
The CO 2 and N deposition-driven land carbon sink is located in northern and tropical forests (Fig. 2). LULCC-induced carbon losses occur across the globe but are most apparent in tropical latitudes and certain hotspot regions across the northern hemisphere (China, USA and West Eurasia), which recently have been (within our study period) densely forested ( Fig. 2 and Supplementary Fig. 1). There are signs of European carbon sink in part driven by a/reforestation or natural land (re)establishment. The impact of changes in climate can also be detected across the globe, with vegetation and soil carbon losses in the Amazon, the Sahel and South Africa, and climate-induced carbon gains in east Brazil, Australia, and across the high northern latitudes ( Fig. 2 and Supplementary Fig. 1). In summary, the DGVM ensemble agrees on the large sink located in the world's forests. However, the models are not in full agreement in the direction of change in ecosystem carbon across much of the globe, in particular in regions with competing CO 2 and LULCC effects (Fig. 2a).

Key processes and uncertainties behind the sink
The top-down global GCB budget constraint combined with the bottom-up DGVM estimates give high confidence the land has been a net sink of carbon over the last 60 years, however, the relative contribution of each driver (e.g. CO 2 , N deposition, climate, LULCC), and ecological attribution (vegetation, soil) of the sink remains elusive.
The DGVM multi-model mean suggests a net global increase in both vegetation (ΔC v = 28 ± 26 (mean ± std. dev) PgC) and soil (ΔC s = 21 ± 32 PgC) carbon stocks, although the uncertainty is large (Fig. 3a, b), with two models (LPX-Bern and YIBs) simulating a net loss of vegetation carbon and three models (CLASSIC-N, ISAM and LPJ-GUESS) simulating a net loss of soil carbon (Fig. 3c).
This relatively high uncertainty is in part due to large opposing fluxes (Fig. 3a) driven by increasing atmospheric CO 2 and N deposition (carbon increase) vs LULCC (carbon decrease at global scale), with soil carbon gains (from rising CO 2 ) also counteracted by negative fluxes driven by climate change (Fig. 3b). In general, there are no observational constraints on long-term (pre-satellite era) changes in global vegetation or soil carbon, and so reducing model uncertainty is challenging. Therefore, although the DGVMs generally capture the global net carbon sink, the attribution to vegetation or soil stocks is completely unconstrained (Fig. 3c) and is driven by alternate model structures and parameterisations between DGVMs.
DGVMs have inherently different baseline productivity, allocation, and turnover rates as well as varying degrees of sensitivity of processes to environmental change. To gain a deeper understanding of the processes driving carbon sink changes, we use our attribution framework (see Methods) to associate the modelled changes in vegetation and soil carbon to changes in inputs (NPP for vegetation and litterfall for soil) and outputs (turnover rates for vegetation and soil). For each of these processes, we identify their changes due to CO 2 , climate change or LULCC, according to the models.
The ensemble-simulated increase in global vegetation carbon over the last 60 years (ΔC v = 28 ± 26 PgC; Fig. 4a) was driven by enhanced plant productivity in response to rising CO 2 concentrations and N deposition, with small (but highly uncertain) losses due to changes in vegetation turnover (ΔĈ v,input,CO2 = 87 ± 26 PgC, ΔĈ v,output,CO2 = −3 ± 21 PgC; Fig. 4a). Although the central estimate of the net CO 2 and nitrogen deposition biomass response appears relatively well constrained at global scale (ΔC v,CO2 = 84 ± 31 PgC), model estimates range from a gain of 30 PgC (YIBs) to a gain of 150 PgC (LPJ) in biomass due to rising atmospheric CO 2 ( Supplementary Fig. 2).
Model spread is driven by a combination of uncertainty in the response of both NPP (e.g. driven by uncertainty in leaf-level photosynthetic response, scaling to canopy and landscape scales, carbon allocation and nutrient limitation) and vegetation turnover (τ v ) (e.g. driven by the treatment of nutrient limitation coupled with differences in allocation and in some cases forest-stand packing constraints 13 ) to rising atmospheric CO 2 and N deposition ( Fig. 5 and Supplementary  Fig. 3). There is some confidence in the direction of change in global and regional NPP, albeit not in the magnitude, whereas models disagree on the direction of Δτ v due to rising CO 2 (Fig. 4). CO 2 and N deposition-driven Δτ v depends on changes in stand dynamics (in some models) and plant allocation. For example, an increase in short-lived root production, at the expense of wood growth, to alleviate nutrient limitations on plant growth can reduce τ v . This increase in root allocation also changes vegetation nitrogen demand (because of the substantially higher C:N in wood compared to roots), which feeds back onto the whole-plant productivity response to rising CO 2 . Conversely, some models increase their wood allocation fraction when NPP (or production) increases, which will increase τ v 14 .
This increase in vegetation stocks has been partially offset by emissions associated with LULCC from vegetation (ΔC v,LULCC = −58 ± 21 PgC; Fig. 4a), predominantly in the tropics (Supplementary Fig. 1), which increases vegetation turnover (loss of biomass to the atmosphere or wood products). About half (ΔĈ v,input,LULCC = −23 ± 22 PgC, 40%; Fig. 4a) of the net LULCC-driven global vegetation carbon losses are from a reduction in inputs (NPP) to the land, as crops or pastures often have lower productivity than the forests that they replace. The other half (ΔĈ v,output,LULCC = −29 ± 18 PgC; Fig. 4a) is attributable to changes in turnover related to forests containing woody biomass, which has a slower turnover than leaves and fine roots. The turnover effect is further aggravated through the phenomenon that due to their longer-lived biomass, forests take up more carbon under rising CO2 levels, but this sink is lost by clearing for agricultural use. This is known as the "loss of additional sink capacity" (LASC) 15,16 and accounts for the impact of environmental (CO 2 and climate) changes on carbon uptake of deforested land (found in S3 simulation) compared to potential vegetation (found in S2 simulation). Ref. 16 attributes~40% of DGVM estimated LULCC vegetation and soil carbon losses to the LASC, albeit with a different methodology and study period.
All models simulate a net loss of global biomass following LULCC, although the magnitude is not well constrained. Model spread is driven by differences in LULCC process representation in models 7 . We find that models that include wood harvest (routine harvest of established managed forests), grazing harvest, or shifting cultivation simulate larger C v losses than models without these land management processes ( Supplementary Fig. 4), in line with earlier studies showing the high importance of land management on vegetation carbon stocks 17 . Further, in northern ecosystems, models do not agree on the direction of ΔNPP or Δτ v following LULCC ( Supplementary Fig. 5a). In addition to the land management processes mentioned above, uncertainty in large-scale regional changes in C v is in part driven by smaller-scale regional compensations, where losses in Russia and USA are somewhat countered by a forest regrowth sink in Eurasia ( Supplementary Fig. 1). Therefore, the balance between these two opposing fluxes determines net northern changes in C v from LULCC, with uncertainties in both carbon loss fluxes and regrowth uptake 18,19 . In addition, spatial patterns of northern LULCC losses and gains are not entirely consistent between models ( Supplementary Fig. 6), which adds to uncertainty in local, regional, and hemispheric net LULCC-driven changes in vegetation carbon.
Changes in climate have a much smaller but similarly uncertain influence on vegetation processes at a global scale, with no agreement among models on the direction of change in biomass driven by shifts in NPP (ΔĈ v,input,CLIM = −1 ± 14 PgC over 1959-2020) or vegetation turnover (ΔĈ v,output,CLIM = 3 ± 9 PgC) at global scale (Fig. 4a). However, models align closer at regional scales, where northern warming has stimulated productivity (ΔĈ v,input,CLIM,North = 11 ± 6 PgC; Supplementary  Fig. 5a), and in the tropics, reductions in productivity (due to warming and/or changes in precipitation) lead to losses of vegetation (ΔĈ v,input,CLIM,Tropics = −12 ± 12 PgC; Supplementary Fig. 5b). Overall, Increased production and loss of biomass enhance soil inputs (ΔĈ s,input,ALL PgC) and drives the growth in global soil carbon stocks (ΔC s PgC), albeit with low confidence in the magnitude (Fig. 4b). Our analysis also indicates that enhanced CO 2 concentrations increase soil turnover, which is a direct result of enhanced inputs to the fastturnover litter and surface soil pools. Our methodology simplifies the TRENDY model structures into a single soil pool to represent the entire soil system, and so a relative increase in the faster pools moves the aggregate soil pool turnover (τ s ) towards that of the litter and surface soil pools-a phenomenon known as false-priming 20 . In addition, changes to soil conditions (e.g. increased soil moisture 13 ) following increased atmospheric CO 2 can impact τ s , although it is difficult to separate the false-priming and actual changes in τ S with the set of simulations used in this study. This simplification leads to low confidence in the partition of net soil carbon changes into the input and turnover-driven changes, whereas the net change is well constrained (ΔC s,CO2 PgC).
Reduced litterfall (input to soil) due to the conversion of forests to agricultural land is the predominant pathway of soil carbon loss, in line with previous studies 21 (ΔĈ s,input,LULCC = −36 ± 76 PgC, ΔC s,LULCC = −25 ± 29 PgC). This is particularly evident in the tropics (Supplementary Fig. 7b) with large-scale deforestation over the last several decades, while deforestation and reforestation are more in balance in the extra tropics. However, there is relatively high uncertainty in the magnitude of LULCC impacts on soil carbon inputs (ΔĈ s,input,LULCC ) (Fig. 4b), with models ranging from losses of over 100 PgC to gains of 45 PgC over the period 1959-2020. Land management processes such as crop and wood harvesting (in particular, the frequency of harvest and fraction of biomass removed and respired elsewhere) are leading uncertainties in simulating LULCC impacts on soil carbon (Supplementary Fig. 4) 22,23 . Turnover-driven changes play less of a role in soil C changes in our models. This could be due to offsetting effects such as a decrease in soil C due to the altered quality of the litter input (a higher fraction of faster decomposing material for agriculture as compared to forests) or a reduction in fire-related losses with the transformation of natural ecosystems 21 .
Increased heterotrophic respiration rates and subsequent reductions in the turnover time of the soil pool due to global warming entirely offset the increase in soil carbon from climate-enhanced NPP and litterfall (Fig. 4b). Overall, changes in climate caused a net loss (ΔC s,CLIM = −13 ± 12 PgC) of global soil carbon. We find the impact of climate change on carbon storage is accelerating, in particular for northern soil and tropical biomass, where carbon stock changes from increased/decreased productivity (for northern soil and tropical  biomass, respectively) have accelerated by 19 TgC yr −2 and 9 TgC yr −2 , since the turn of this century (Supplementary Fig. 8). This equates to annual mean carbon gain/losses due to changes in productivity of 1.6 PgC yr −1 and −0.3 PgC yr −1 over 2001-2020 compared to 0.6 PgC yr −1 and −0.2 PgC yr −1 over 1961-2000, for northern soil and tropical biomass, respectively ( Supplementary Fig. 8b, e). This acceleration is likely driven by warming in the north (2001-2020 April-August mean 0.9°C warmer than 1961-2000) and the tropics (2001-2020 annual mean 0.6°C warmer than 1961-2000) 24,25 , as well as potential changes in sensitivity of the biosphere to warming 26 .

Discussion
The 'observed' GCB global net land sink is captured by the TRENDYv10 DGVM ensemble. We find a robust agreement that rising atmospheric CO 2 concentrations and N deposition drive the land carbon sink, and changes in climate along with LULCC lead to a source of carbon in the atmosphere. Further, the DGVMs corroborate that the net land sink is located in northern latitudes, with the tropics more carbon neutral, in line with recent independent top-down estimates 3,8 .
However, the models do not entirely agree on the partition of the net land sink between vegetation and soil. Thirteen of 18 models indicate an increase in global vegetation and soil carbon stocks, although the magnitude of change is highly unconstrained (ranging from~0 to +80 PgC over the past 60 years). The remaining five models simulate a reduction in either vegetation or soil stocks, highlighting the disparity between models regarding internal carbon cycling processes.
There are no observations of global carbon stocks covering the entire study period, but a recent synthesis of observation-based biomass change estimates in global forests suggests a net biomass sink of 0.3-2.1 PgC yr −1 since 2000 27 , giving confidence to the models that simulate net gains. However, these observational studies do not constrain the magnitude of the vegetation sink any more than the DGVM mean biomass sink of −0.1 to 1.7 PgC yr −1 (averaged over 2001-2020). Furthermore, we cannot even be certain in the direction of change in global soil carbon as direct observations of large-scale changes in soil carbon stocks do not exist 2,4 .
Our analysis indicates that baseline vegetation turnover rates and the turnover response to rising atmospheric CO 2 and N deposition are the key uncertainties in modelled ΔC v , with a smaller but still significant contribution from changes in NPP. Differences in modelled vegetation turnover are driven by variations in whole-plant mortality rates, simulated tissue lifespan, and the allocation of NPP to plant components with inherently different turnover rates and their responses to changing environmental conditions 14,28 . For example, some models may reduce plant mortality and turnover as atmospheric CO 2 concentrations rise. Increased carbon stores could supply maintenance respiration in periods of photosynthetic stress 29,30 , and increases in water-use efficiency may alleviate the impact of drought 31,32 . In contrast, shifts in plant community composition or selfthinning dynamics may increase stand-level mortality 33 . These structural differences lead to no agreement on the direction of Δτ v in the models presented here.
In addition to turnover-related uncertainties, the increase in NPP and biomass due to rising CO 2 is not well constrained. The choice of photosynthesis and stomatal model determines the leaf-level response, with a wide range of implementations and outcomes in DGVMs 34 . Scaling from leaf to canopy level causes further discrepancies due to different structural assumptions between models as to the vertical distribution of light, nitrogen, photosynthetic capacity, and the treatment of sunlit and shaded leaves 34 . Further, NPP and biomass production is constrained by stoichiometric nutrient requirements. Nitrogen availability mediates the biomass response to rising CO 2 35 , and therefore it is likely that differences in the inclusion and representation of carbon-nitrogen coupling between models contribute to uncertainty in ΔC v 36,37 . Overall, these structural uncertainties manifest as a four-fold range in estimates of the CO 2 effect on NPP and biomass (gains of 41 PgC (CLASSIC-N) up to 157 PgC (ISAM) over the previous 60 years).
Furthermore, LULCC-driven reductions in NPP and woody carbon stocks have removed a substantial portion of global vegetation carbon. However, the exact magnitude of carbon loss is not well constrained by the DGVMs. Biases in the carbon density of deforested land 38 and model-specific choices of whether to replace forest or grassland for agricultural expansion are leading sources of error when calculating LULCC losses 39 . Our results highlight that, firstly, variation in the representation of relevant processes, in particular, forest management (e.g. treatment of wood harvest) and shifting cultivation, leads to systematic differences between modelled ΔC v 7,40 . Second, DGVMs are missing carbon losses associated with degraded forests, which can exceed deforestation losses, indicating simulated LULCC losses could be underestimated [40][41][42] . Finally, uncertainty exists with LULCC maps used to drive DGVMs as historical land-use is not perfectly known and products differ in the land cover types and transitions included 40 .
The models suggest that soil has been a net carbon sink globally over the past six decades, although there is over 100% uncertainty on the magnitude (21 ± 32 PgC). This relatively large range is due to the opposing impacts of rising CO 2 and climate/LULCC effects which partly cancel out. Further, the impact of LULCC on soil carbon is also difficult to quantify robustly with the suite of DGVMs used here, as models do not agree on the sign of Δf vs or Δτ s following LULCC. In general, changes in soil carbon stocks are difficult to constrain due to a lack of large-scale observations. It is likely historical LULCC caused additional carbon fluxes as a result of land erosion 43,44 , degradation of agricultural soils 45,46 , and losses from drained peatlands 47 . These management processes are not included in models and most do not simulate peatlands at all, potentially causing errors in simulated soil carbon fluxes.
Overall, the soil carbon sink is caused by increased litter inputs from CO 2 -driven vegetation growth. There is some evidence additional litter and root exudates enhance soil carbon stocks when nitrogen availability is high 48 . It is important to stress that the model simulations include widespread nitrogen deposition, which may have helped to sustain a strong CO2 response of biomass and soil carbon sinks 5,49 .
However, additional soil inputs can accelerate organic matter decomposition to release plant-available nitrogen (via priming effects 50 ), which reduces soil carbon storage [51][52][53] . DGVMs do not simulate priming effects or explicitly account for microbial activity, mineral association, aggregation, or mycorrhizal fungi interactions, all of which regulate soil carbon turnover and accumulation rates 54 . DGVMs represent soil carbon with a set of cascading pools, where decomposition losses are determined by first-order decay rates (although this is an active area of model development 55 ). Therefore, the modelled soil sink has to be interpreted in the context of model structure and the limited ability of DGVMs to capture governing soil processes.
In general, it is difficult to attribute processes (soil inputs vs turnover) to changes in soil carbon for two main reasons. First, soil model structure leads to false-priming effects 20 . The apparent reduction in soil turnover is a consequence of increased surface soil carbon relative to deep soil, as opposed to a change in actual turnover rates. Second, our process attribution methodology treats the soil as a single pool. Therefore, the calculation of τ s ( C s R h ) may include deep, inactive (on decadal timescales) soil carbon, and underestimate the turnover rate of the active soil. Hence, our process attribution of soil carbon changes is not well constrained. For example, we estimate a 300 PgC increase in soil carbon from increased inputs for CLM5.0, a result which is in part driven by the large carbon stocks in high northern latitudes leading to a large estimate of τ s 56 . We do account for this bias by scaling output to actual changes in modelled carbon stocks (see Article https://doi.org/10.1038/s41467-022-32416-8 methodology), but our results highlight the urgent need for evaluation of individual carbon pools and fluxes (leaf, root, wood and each soil pool) separately, and not bulk terms.
In addition to the CO 2 -driven sink, evidence exists for a substantial impact of forest regrowth on the northern sink over the previous 60 years 18,19 . DGVMs have some capacity to capture regrowth effects (predominantly following land abandonment), however, most do not simulate forest demography or detailed forest management (e.g. enhanced stock densities 57 ) and so could underestimate the actual regrowth carbon sink 58 . Further, not all DGVMs simulate disturbance in unmanaged land (e.g. wind and pests), and therefore, natural disturbances beyond fire and the subsequent regrowth will not be captured by the DGVMs here 18 . Moreover, not all models simulate fire mortality, and fire models do not always capture observed spatiotemporal patterns and associated carbon emissions 59 , as well as large uncertainties due to issues with forcing data 60 .
In addition to direct human-caused carbon losses, alternative drivers of mortality can reduce carbon stocks. For example, climateinduced mortality, combined with a shorter tree lifespan resulting from faster growth 61,62 , currently weaken the Amazon forest sink 63 . However, DGVMs show no sign of increased carbon turnover due to changes in climate or faster growth, although this discrepancy is not unexpected as DGVMs do not include detailed mortality processes, e.g. drought-mortality 64 , although this is an area of current model development 65,66 . The DGVMs do indicate a decline in tropical productivity due to climate change, in line with remote-sensing and upscaled in-situ observations 67 , indicating current temperatures may exceed tropical plant thresholds 24 .
Here we have shown that the long-term evolution of the global net land sink is well estimated by an ensemble of state-of-the-art DGVMs. However, a complete process-based understanding is still lacking due to several model shortcomings. Specifically, above/below-ground partitioning of the sink is highly uncertain due to too simplistic representation of internal carbon cycling. Model improvement into plant allocation, tissue lifespan and mortality, as well as the inclusion of process-based soil carbon and nutrient cycling, should be a high priority moving forward.

TRENDYv10 models
We analyse output from 18 DGVMs that are part of a recent model intercomparison project, TRENDYv10 and the Global Carbon Budget 2021 (GCB2021) 8  The models are forced with a merged monthly Climate Research Unit (CRU) 70 and 6-hourly Japanese 55-year Reanalysis (JRA-55) 71 data set. The models are also forced with atmospheric CO 2 72 , gridded nitrogen deposition 73 and nitrogen fertiliser 74 . DGVMs use the HYDE (v3.3) land-use change data set 75,76 , which provides annual pasture and cropland areas at a global scale, and includes improvements in the spatial distribution of agricultural regions 77 . Several models (CABLE-POP, CLM5.0, JSBACH, LPJ-GUESS, LPJ, and VISIT) also use harmonised land-use change data (LUH2-GCB2021 8 ), which provides information on sub-grid-scale land-use transitions.
To isolate the response of the land to each driver (CO 2 , climate, LULCC), each model performs four simulations: S0 (fixed preindustrial atmospheric CO 2 and land-use, recycled 1901-1920 climate), S1 (transient atmospheric CO 2 , recycled 1901-1920 climate, and fixed pre-industrial land-use), S2 (transient atmospheric CO 2 , transient climate, and fixed pre-industrial land-use), and S3 (transient atmospheric CO 2 , transient climate, and transient industrial land-use). Nitrogen deposition varies temporally in simulations S1-S3. Therefore, the transient CO 2 (+N deposition) effect on the 'natural' land sink is calculated by S1−S0, S2−S1+S0 is the climate effect on the 'natural' land sink, S3−S2 is the LULCC effect, and S3 is the net effect. There exists an artefact in the HYDE3.3 data causing a large land-use transition and emission peak around 1960. To correct this, we replace the LULCC estimates for 1959-1961 with the average of 1958 and 1962 in each DGVM.
Net Biome Productivity (NBP) from the S3 simulation represents the net land sink and can therefore be compared to the 'observed' Global Carbon Budget land constraint, which is calculated in the GCB as fossil fuel emissions-atmospheric carbon growth rate-ocean carbon sink (see Fig. 1 in the main text and Table 5 in ref. 8 and data taken from https://doi.org/10.18160/gcp-2021). Atmospheric CO 2 concentration measurements began in 1959 78 , and so this independent constraint on the land sink covers the period 1959-2020, which defines the study period used throughout our analysis.

Data processing
First, we calculate for each model the global and regional (two regions: north of 30°N and south of 30°N) annual mean values for NBP, net primary productivity (NPP), heterotrophic respiration (R h ), vegetation carbon (C v ) and soil carbon (C s ), which includes litter (C litter ) and coarse woody debris (C CWD ) for 1959-2020. We exclude the product pool from our analysis as we do not have carbon emission data from wood products available.

Process attribution framework
To attribute changes in C pools to different processes, we quantify C pool changes attributable to changes in productivity and changes in turnover time. We do this using an analytical approximation of C pool dynamics corrected for non-steady state and more complicated behaviour of the models. First, we approximate the steady-state C pool size given the inputs and turnover time of a given pool in a given year. We then quantify the change in the steady-state pool sizes estimated for 1959 and the year in question and partition this change proportionally among inputs and turnover time. Finally, to adjust to nonsteady state conditions and more complicated model structures (i.e. multiple pool dynamics, not a single pool as in this approximation), we adjust the attributed proportions by adding the difference in actual modelled change in carbon stock and the steady-state approximation, scaled by the relative magnitude of productivity and turnover-driven changes.
The change in carbon pools can be expressed as: where we define dC dt as the change in carbon pool per year. To quantify the internal processes impacting the output flux, we can define OUTPUT = C τ , where τ is the mean turnover time (the inverse of the mean turnover rate) of the carbon pool, C. Therefore, we can rewrite Eq. 1 as: We separate the land into two main pools, vegetation (C v ) and soil (C s ). C v is a direct model output and includes leaf, root and wood biomass. Similarly, C s ,C litter ,C CWD are direct model outputs, which we sum together for our analysis. The rest of the methodology, C s refers to the sum of these three pools. Not all models simulate litter and CWD pools. ISAM does not simulate litter and only CABLE-POP, CLM5.0, DLEM and LPJ-GUESS simulate CWD. Further, for CLM5.0, the litter pool is incorporated in the soil pool output, and so is also excluded from the sum to avoid double counting.
We can now specify Eq. 2 for vegetation and soil as: where f vs is the carbon flux from vegetation to soil, from either litterfall, mortality or direct transfer from roots to soil. We calculate annual τ v (using Eq. 3), τ s and f vs using model output as follows: where ΔC v and ΔC s are the annual changes in vegetation and soil carbon. We want to be able to define ΔC in terms of input fluxes and turnover times only, as these are the two aggregated mechanisms driving changes in the pools. To do this, we use the fact that the pools and the input and output fluxes are generally much larger than the change in carbon pools each year to make the assumption that ΔC ' 0. Rearranging Eqs. 3 and 4, we can estimate the steady-state carbon pools in any particular year as: whereĈ v andĈ s are steady-state approximations of the actual carbon pools C v and C s . We can now define the change in carbon pools over our 60-year study period, 1959-2020, in terms of inputs and turnover time: ΔC v ' ΔĈ v = NPP 2020 τ v,2020 À NPP 1959 τ v,1959 ð7Þ ΔC s ' ΔĈ s = f vs,2020 τ s,2020 À f vs,1959 τ s,1959 : In order to separate the impacts of changes in inputs and turnover times, we next define NPP 2020 = NPP 1959 + ΔNPP, f vs,2020 = f vs,1959 + Δf vs , and τ 2020 = τ 1959 + Δτ. We now substitute these into Eqs. 7 and 8: The right-hand side of the equation contains three terms corresponding to the change in carbon storage due to changes in inputs (NPP for vegetation and litterfall for soil), outputs (bulk turnover), and an interaction term. Our approach, therefore, extends the factor separation approach by ref. 21, which applied it to attribute simulated soil carbon changes into input-driven, turnover-driven change and a synergy term to also cover vegetation carbon changes.
The change in vegetation carbon due to each term can be written as: ΔĈ v,input = ΔNPPτ v,1959 , ΔĈ v,output = NPP 1959 Δτ v , ΔĈ v,interaction = ΔNPPΔτ and ΔĈ v = ΔĈ v,input + ΔĈ v,output + ΔĈ v,interaction , with equivalent terms for soil carbon.
Finally, we adjust the one-pool steady-state estimates for the difference in the model simulations that arise due to the combination of steady-state assumption (Eqs. [5][6][7][8], the grouping of all pools into single vegetation and soil pools (as the distribution between component pools changes over time), and the linearisation not capturing the full dynamics of land carbon cycling (Eqs. 9 and 10). Therefore, the estimates of changes in carbon pools (ΔĈ v and ΔĈ s ) do not equal the actual change in C pools (ΔC v and ΔC s ), which are a direct model output (see Supplementary Fig. 9).
To address this issue, we calculate for each model the difference between the actual change in carbon stocks and our approximation (δ v = ΔC v À ΔĈ v and δ s = ΔC s À ΔĈ s ). We make the simple assumption that productivity and turnover are proportionally responsible for the mismatch between steady state and the actual carbon stock, and so adjust each term on the right-hand side of Eqs. 9 and 10 by the relative magnitude of each term multiplied by δ v or δ s for vegetation and soil, respectively. For example, ΔĈ v,input is adjusted by adding: We calculate Eqs. 9 and 10 for each of the eighteen models, apply the appropriate scaling and use them to produce Fig. 4.

Uncertainty analysis
The spread of model estimates of changes in carbon stocks and fluxes is used to determine the uncertainty of our results. For land flux estimates (NBP) and ΔC v and ΔC s , we show 'mean ± std. dev.'. We decompose the uncertainty (defined as the model spread) in ΔĈ v and ΔĈ s into four components: baseline input (NPP 1959 or f vs,1959 ), change in inputs (ΔNPP or Δf vs ), baseline turnover (τ v,1959 or τ s,1959 ), and change in turnover (Δτ v or Δτ s ). For each of these terms in Eqs. 9 and 10, we calculate the standard deviation in ΔĈ v and ΔĈ s using the multimodel mean values of all other terms in the equations and the individual model values for that term.

Data availability
All code and post-processed data generated in this study are available at https://doi.org/10.5281/zenodo.6884342. The raw model output is available at the following sftp site: trendy-v10@trendy.ex.ac.uk. Access will be granted by first contacting Stephen Sitch (S.A.Sitch@exeter.ac.uk).